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1.  TRANSIENT  FLOW  SIMULATION 


The  U.S.  Army  Research  Laboratory  (ARL)  uses  a  variety  of  computational  fluid  dynam¬ 
ics  (CFD)  techniques  for  the  numerical  simulation  of  air  blast. ^  High  resolution  simulations  of 
blast  flows  are  accomplished  through  the  use  of  two-dimensional  (2-D)  and  three-dimensional 
(3-D)  Euler^  and  Navier-Stokes^  flow  solvers.  Typical  2-D  flow  simulations  require  tens  of 
hours  of  computer  time  to  complete,  while  3-D  calculations  can  require  hundreds  of  hours.^ 
With  high  resolution,  multi-dimensional  simulation  come  huge  data  storage  requirements 
and  time-consuming  data  analysis.  Some  blast  studies  require  the  capability  to  estimate 
trends  in  flow  field  characteristics  through  parametric  analysis.  Often,  this  type  of  analysis 
consists  of  a  large  number  of  calculations  that  must  be  completed  in  a  short  period  of  time 
to  address  critical  project  issues.®  To  address  this  need,  a  flow  solver  has  been  developed  to 
simulate  transient,  compressible  fluid  flow  in  a  one-dimensional  (1-D)  computational  domain 
with  variable  area.  This  new  CFD  tool  is  called  the  Transient  One-Dimensional  Flow  Solver 
(TOFS).  A  description  of  the  TOFS  algorithm,  code  features,  and  the  presentation  of  some 
validation  test  cases  form  the  basis  of  this  report. 


2.  GOVERNING  EQUATIONS 


The  motion  of  a  fluid  can  be  completely  described  by  the  conservation  of  mass,  momen¬ 
tum  and  energy.  Conservation  of  momentum  for  an  inviscid,  Newtonian  fluid  is  described 
by  the  Euler  equation.®  When  gravitational  effects  are  ignored,  the  Euler  equation  takes  the 


general  form 


Dui  _  ^ 

^  Dt  dxi 


(1) 


in  which  Xi  represents  the  spatial  dimensions,  and  Ui  represents  the  velocity  components  in 
these  dimensions. 


When  applied  in  only  one  dimension,  the  conservation  laws  take  the  form  of  the  partial 
differential  equation^  (PDE) 

dt  dx 


(2) 


in  which 


q  = 


pA 

pAu 

pA{e  +  _ 


(3) 


and 


F  = 


pAu 

pAv?  -I-  pA 
u[pA{e  +  +p] 


(4) 


1 


In  this  form,  p  is  the  fluid  density,  A  is  the  area,  u  is  the  velocity,  p  is  the  pressure, 
and  e  is  the  specific  internal  energy.  For  problems  of  interest  to  blast  modeling,  the  fluid  is 
considered  to  be  an  ideal  gas,  for  which  the  specific  internal  energy^  is  given  by  Equation  5. 


_  P 
®  p(7  -  1) 


(5) 


In  this  relationship,  7  represents  the  ratio  of  specific  heats  and  is  assumed  to  be  constant 
for  the  range  of  temperatures  that  are  expected  in  shock  tube  applications. 


3.  SOLUTION  METHOD 


Time-accurate  simulation  of  blast  flow  in  air  is  accomplished  here  by  first  transforming 
the  problem  geometry  into  a  computational  domain  where  grid  points  may  be  uniformly  or 
variably  spaced  along  the  length  of  the  domain.  The  governing  equations  are  discretized 
and  then  solved  iteratively  to  march  the  solution  through  time.  The  MacCormack  finite 
difference  scheme  is  used  as  the  solution  method.  This  explicit  method  uses  a  predictor- 
corrector  approach  for  advancing  the  solution  from  the  current  time  step  to  the  next.  In  the 
predictor  step,  a  forward  difference  is  used  to  compute  an  interim  solution  at  each  grid  point 
as  illustrated  in  Equation  6. 


Here,  the  qij  terms  denote  the  conservation  parameter  from  the  vector  of  dependent  vari¬ 
ables,  q,  from  Equation  3,  in  which  i  represents  the  specific  conservation  parameter  of  mass, 
momentum  and  energy,  and  j  represents  the  point  in  the  grid  at  which  the  solution  is  being 
performed.  Likewise,  the  Fij  and  Fij+i  terms  are  used  to  identify  the  particular  independent 
variable,  i,  at  grid  points  j  and  j  -h  1,  respectively.  The  asterisk  (*)  is  used  to  identify  the 
predicted  conservation  parameter  and  n  is  used  to  identify  variables  taken  at  the  current 
point  in  time.  The  value  of  the  time  step  is  given  by  At,  and  Aajj  is  the  inter-grid  spacing 
at  grid  point  j. 

After  the  predictor  step  is  used  to  compute  interim  values  of  the  q  vector  at  every  grid 
point,  the  corrector  step  uses  the  predicted  values  and  a  backward  difference  to  compute  the 
solution  at  the  next  time  step  for  each  grid  point,  as  shown  in  Equation  7. 


To  account  for  area  changes  in  the  computational  domain,  the  momentum  equation,  the 
second  elements  of  Equations  3  and  4,  is  predicted  in  the  manner  shown  in  Equation  8. 


At 

Axj 


(8) 
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And  the  corrector  step  takes  the  form 


jpn 

•^y-1 


p 


(9) 


After  the  vector  of  dependent  variables  has  been  updated  and  the  solution  advanced  in 
time,  a  new  time  step  must  be  computed.  The  values  of  the  q  vector  are  used  to  compute 
the  sound  speed  and  the  particle  velocity  at  each  grid  point.  The  inter-grid  spacing  is 
then  divided  by  the  sum  of  the  sound  speed  and  the  particle  velocity  to  compute  a  signal 
propagation  time  for  each  grid  point.  The  minimum  signal  propagation  time  in  the  grid  is 
multiplied  by  a  stability  factor^  which  is  less  than  1.0  and  the  resulting  value  becomes  the 
next  time  step.  This  method  of  computing  the  time  step  is  summarized  in  Equation  10 

At  =  CFL  *  MINj  ( 7-4^)  (19) 

\\Uj\  +  CjJ 

in  which  CFL  is  the  stability  factor,  Axj  is  the  grid  spacing  at  grid  point  j,  is  the 
absolute  value  of  the  velocity  at  grid  point  j,  and  Cj  is  the  sound  speed  at  grid  point  j. 


4.  NORMAL  SHOCK  WAVES 

Transient  flow  in  shock  tubes  is  characterized  by  the  propagation  of  a  shock  front  through 
undisturbed  air.  This  leading  shock  may  be  followed  by  transient  flow  dictated  by  the 
geometry  of  the  shock  tube  or  other  secondary  shocks.  In  most  shock  tubes,  the  primary 
direction  of  flow  is  parallel  to  some  reference  axis.  When  the  leading  shock  is  perpendicular 
to  the  flow  direction,  it  is  referred  to  as  a  “normal  shock  wave.”  If  the  conditions  of  the 
undisturbed  gas  are  known,  then  the  conditions  of  the  fluid  immediately  behind  the  normal 
shock  can  be  completely  described  if  only  one  parameter  of  the  shock  wave  is  known.  This 
information  can  be  used  to  compare  a  variety  of  fluid  parameters  from  blast  experiments  to 
those  values  generated  by  a  computational  prediction,  and  this  is  useful  in  code  validation 
efforts.  The  ^nkine-Hugoniot  Equations  relate  the  gas  parameters  across  the  shock  wave, 
and  are  derived  from  the  conservation  of  kinetic  energy  between  the  states  on  either  side  of 
the  shock.  ^9’  ^ 

The  speed  of  the  shock  wave,  the  particle  velocity  behind  the  shock,  and  the  density 
behind  the  shock  are  related  to  the  pressure  ratio  across  the  shock  by  the  relationships  in 
Equations  11,  12,  and  13.  Here,  p^,,  Poo,  and  Coo  are  the  ambient  atmospheric  pressure, 
density,  and  sound  speed,  respectively.  The  shock  speed  is  represented  by  the  symbol  Ug, 
while  the  pressure,  density,  and  velocity  behind  the  shock  are  identified  as  pi,  pi,  and  «i, 
respectively. 


^  =  1. 
cL  27 


(7-1)  +  (t  +  1)^ 


(11) 
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Coo 


7  -  1  Pi  I  7  +  1 
7  VPoo  )[  27  Poo  27  _ 

^  (7-  1)  Poo  +  (7  +  1)  Pi 
Poo  (7+ l)Poo  +  (7-  l)Pl 


(12) 

(13) 


5.  BLAST  FLOW  MEASUREMENTS 


Before  presenting  the  comparison  of  shock  tube  experimental  data  with  computational 
results,  it  is  necessary  to  discuss  the  types  of  flow  data  that  can  be  obtained  in  such  ex¬ 
periments  and  how  this  information  can  be  used  to  derive  other  fluid  parameters.  The 
primary  means  of  collecting  flow  history  data  in  shock  tubes  is  through  the  use  of  pressure 
measurements.  By  measuring  the  static  and  stagnation  overpressures  at  a  point  in  the  flow 
and  knowing  the  ambient  atmospheric  conditions,  it  is  possible  to  determine  some  useful 
information  about  the  blast  event.  The  dynamic  pressure  of  a  fluid  in  motion  represents  the 
kinetic  energy  per  unit  volume  of  the  fluid  and  is  a  function  of  the  fluid  density  and  velocity. 
When  the  fluid  being  considered  is  an  ideal  gas,  the  dynamic  pressure  can  also  be  expressed 
as  a  function  of  the  pressure  and  the  local  Mach  number  which  may  be  obtained  from  the 
Pitot  Equations.  For  subsonic  flow,  the  dynamic  pressure,  q,  may  be  calculated  from  the 
static  pressure,  p,  and  the  stagnation  pressure,  po,  according  to  Equation  14. 


Q 


7  2p 

27-1 


(14) 


For  low  speed  subsonic  blast  flows,  small  errors  or  oscillations  in  the  measured  static 
and  stagnation  overpressure  histories  can  lead  to  difiiculty  in  accurately  calculating  the  local 
Mach  number  and  dynamic  pressure  behind  these  weak  shocks.  To  reduce  the  error  in 
calculating  dynamic  pressure  in  such  flows,  a  differential  pressure  gauge^  can  be  used.  This 
device  has  a  single  pressure-sensing  element  that  directly  measures  the  difference  between 
stagnation  and  static  pressure.  When  a  static  gauge  and  a  differential  pressure  gauge  are 
used  to  measure  the  flow  histories  at  a  common  point  in  the  flow  field,  the  two  records  may 
be  summed  to  yield  the  stagnation  pressure  that  can  then  be  used  to  determine  the  local 
Mach  number  and  dynamic  pressure. 

For  cases  in  which  the  differential  pressure  gauge  and  the  static  pressure  gauge  are  not 
close  enough  to  justify  summation  to  determine  the  stagnation  overpressure,  the  measured 
differential  pressure  may  be  used  to  approximate  the  dynamic  pressure.  Figure  1  illustrates 
the  relationship  of  difierential  pressure  and  dynamic  pressure  to  local  Mach  number.  For 
all  the  cases  presented  in  this  report,  the  local  Mach  number  is  less  than  0.7.  At  this  Mach 
number,  the  ratio  of  differential  pressure  to  ambient  static  pressure  is  0.38  and  the  ratio 
of  dynamic  pressure  to  ambient  static  pressure  is  0.34.  The  relative  difference  between  the 
differential  pressure  and  dynamic  pressure  at  this  Mach  number  is  12%.  This  small  relative 
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difference  is  evidence  that  the  differential  pressure  can  be  reliably  used  to  approximate  the 
dynamic  pressure  for  subsonic  blast  flows. 


Figure  1.  Differential  and  Dynamic  Pressure  as  a  Function  of  Shock  Strength. 

Having  discussed  the  relationships  between  experimental  flow  measurements,  it  is  neces¬ 
sary  to  develop  a  strategy  which  will  be  used  in  employing  experimental  data  to  derive  other 
flow  parameters  and  a  set  of  rules  for  comparing  experimental  data  to  computed  results.  The 
following  is  the  list  of  rules  that  are  used  to  determine  which  flow  parameters  to  compare 
with  computational  results: 

1.  If  only  one  pressure  gauge  is  used  at  a  given  location,  its  results  will  be  directly  compared 
to  the  simulation. 

2.  If  any  two  different  types  of  pressure  gauges  are  used  at  a  common  location,  the  data 
recorded  by  those  gauges  will  be  used  to  derive  the  dynamic  pressure  and  Mach  number 
histories  at  that  location.  If  static  overpressure  is  not  one  of  the  measured  histories,  it 
will  be  derived  from  the  measured  data  as  well. 

3.  For  locations  that  have  at  least  two  different  types  of  gauges,  the  static  overpressure, 
dynamic  pressure,  and  Mach  number  histories  will  be  compared  to  the  simulation. 

4.  For  locations  that  have  only  one  type  of  gauge,  that  information  will  be  compared  with 
the  simulation. 
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6.  SINGLE  STAGE  SHOCK  TUBE  FLOW  SIMULATIONS 


The  intended  use  of  the  flow  solver  described  in  this  report  is  the  numerical  simulation 
of  time-dependent  flow  in  shock  tubes  and  blast  simulators.  This  type  of  analytical  tool  can 
be  used  to  study  possible  shock  tube  configurations  so  that  particular  flow  patterns  can  be 
produced,  experiments  can  be  designed  for  blast  simulators,  predicted  environments  can  be 
predicted  for  gauge  ranging  in  blast  tests,  or  1-D  blast  wave  histories  can  be  generated  for 
use  as  an  inlet  boundary  condition  for  a  separate  numerical  simulation. 

To  validate  the  code  for  modeling  the  types  of  flow  fields  encountered  in  shock  tube 
experiments,  shock  tubes  of  various  configurations  were  selected  and  calculations  were  per¬ 
formed  to  match  experiments  from  each  facility.  These  validation  cases  are  presented  in 
order  of  increasing  geometric  complexity. 


6.1.  Constcint  Area  Shock  Tube 

The  simplest  shock  tube  geometry  is  that  of  a  straight  tunnel  with  constant  cross- 
sectional  area.  One  such  facility  is  the  57-cm  shock  tube  at  ARL.i2This  shock  tube  consists 
of  a  cylindrical  driver  tube  and  expansion  section,  a  square  section,  and  an  additional  cylin¬ 
drical  expansion  region  down  stream  from  the  test  section.  The  cylindrical  tube  sections  are 
57.47  cm  in  diameter,  and  the  square  test  section  is  50.80  cm  on  a  side  to  yield  a  constant 
cross-sectional  area  of  0.26  m^.  The  length  of  the  driver  section  is  adjustable  from  30.48  cm 
to  16.76  m.  The  distance  from  the  beginning  of  the  expansion  tunnel  to  the  measurement 
station  is  30.48  m. 

A  calculation  was  performed  to  model  a  test  of  this  facility  in  which  the  driver  length  was 
91.44  cm,  the  driver  overpressure  was  379.2  kPa,  The  ambient  pressure  was  102.1  kPa,  and 
the  ambient  temperature  was  294.5  K.  The  driver  gas  was  supplied  by  air  compressors  and  the 
pressurized  air  was  at  the  same  temperature  as  the  surrounding  atmosphere.  The  calculation 
was  performed  using  a  computational  domain  which  was  discretized  with  a  uniform  grid 
spacing  of  4  cm.  The  static  overpressure,  dynamic  pressure,  and  Mach  number  histories 
produced  by  the  simulation  are  compared  to  the  measured  data  in  Figures  2  through  4,  and 
Table  1  provides  a  summary  the  comparison  between  the  experiment  and  the  calculation. 

The  static  overpressure  comparison  shows  some  oscillation  in  the  computed  history  im¬ 
mediately  after  the  shock  arrives  at  the  measurement  position.  After  this  period  of  oscilla¬ 
tion,  which  lasts  about  2  ms,  the  computed  history  follows  the  experimental  data  to  within 
2  kPa.  These  characteristics  of  initial  oscillation  after  the  shock  arrival,  followed  by  close 
agreement  with  the  experimental  data  are  evident  in  the  dynamic  pressure  and  Mach  number 
histories  as  well. 

The  oscillations  occurring  near  the  shock  arrival  are  the  result  of  dispersion-the  nu¬ 
merical  propagation  of  plane  waves  through  the  computational  domain.  The  presence  of 
these  oscillations  makes  it  difiicult  to  determine  the  amplitude  of  the  incident  shock  in  time- 
dependent  flow  histories  such  as  those  shown  in  Figures  2  and  3.  However,  the  pressure 
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Dynamic  Pressure  (kPa)  Static  Overpressure  (kPa) 


Figure  2.  57-cin  Shock  Tube  Test:  Static  Overpressure. 


Figure  3.  57-cm  Shock  Tube  Test:  Dynamic  Pressure. 
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Figure  4.  57-cm  Shock  Tube  Test:  Mach  Number. 

immediately  behind  the  shock  can  be  determined  from  the  computed  speed  of  the  shock 
wave.  By  placing  additional  monitoring  stations  in  the  computational  domain  in  the  vicinity 
of  the  target  location,  the  shock  wave  speed  may  be  determined  by  observing  the  difference 
in  shock  arrival  times  at  the  different  positions.  For  this  calculation,  the  shock  wave  speed 
was  found  to  be  436.7  m/s.  Using  the  Rankine-Hugoniot  equations  described  earlier,  the 
incident  shock  overpressure  computed  by  the  code  is  found  to  be  71.2  kPa,  7.4%  greater 
than  the  measured  overpressure.  The  dynamic  pressure  and  local  Mach  number  behind  the 
incident  shock  were  similarly  over-predicted  by  the  computation. 


Table  1.  Summary  of  Results  for  57-cm  Shock  Tube  Test. 


Experiment 

Simulation 

Error 

Shock  Wave  Speed 

Incident  Static  Overpressure 

Dynamic  Pressure  Behind  Shock 

Mach  Number  Behind  Shock 

66.3  kPa 
13.6  kPa 
0.34 

464  m/s 
71.2  kPa 
16.6  kPa 
0.37 

-^7.4% 

+22.1% 

+8.8% 

Static  Overpressure  Impulse  at  100  ms 
Dynamic  Pressure  Impulse  at  100  ms 

0.88  kPa-s 
0.122  kPa-s 

0.90  kPa-s 
0.125  kPa-s 

+2.3% 

+2.5% 
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6.2.  Variable  Area  Shock  Tube 


The  next  level  of  geometric  complexity  to  be  studied  in  the  validation  of  the  code  is  that 
of  a  shock  tube  with  a  varying  cross-sectional  area.  One  such  facility  is  the  ARL  2.44-m 
blast  simulator, a  schematic  diagram  of  which  is  presented  in  Figure  5.  This  facility  has  a 
cylindrical  driver  91.44  cm  in  diameter  and  9.42  m  in  length.  A  converging  nozzle,  located 
at  the  downstream  end  of  the  driver,  channels  the  flow  into  a  throat  section  which  is  half  the 
diameter  of  the  driver  section.  Two  diaphragms  om  the  throat  section  initially  separate  the 


All  Dimensions  in  Meters 
Not  to  Scale 

Figure  5.  2.44-m  Blast  Simulator  Schematic  Diagram. 

high  pressure  driver  gas  from  the  ambient  air.  Before  a  test,  the  pressure  is  increased  from 
ambient  to  the  high  pressure  driver  condition  by  incremental  pressure  increases  across  the 
diaphragms.  The  net  pressure  difference  from  the  driver  condition  to  the  ambient  state  is 
greater  than  either  of  the  diaphragms  can  maintain,  but  the  pressure  differential  across  either 
of  the  diaphragms  is  less  than  the  rupture  pressure  of  the  diaphragms.  Flow  is  initiated  by 
releasing  the  gas  between  the  two  diaphragms,  causing  the  upstream  diaphragm  to  rupture 
first,  followed  by  the  downstream  diaphragm.  At  the  downstream  end  of  the  throat,  there  is 
an  immediate  transition  to  the  expansion  tunnel,  which  is  2.44  m  in  diameter  and  71.78  m  in 
length.  Static  overpressure  and  differential  pressure  histories  are  recorded  in  the  expansion 
section  at  a  position  15.86  m  down  stream  from  the  beginning  of  the  expansion  section. 

There  are  two  significant  differences  between  the  computational  model  of  the  2.44-m 
blast  simulator  and  the  actual  facility.  As  stated  earlier,  the  momentum  calculation  in  the 
MacCormack  predictor-corrector  technique  contains  an  additional  term  to  account  for  the 
rate  of  area  change  per  unit  length  in  the  computational  domain.  In  the  case  of  the  2.44-m 
blast  simulator,  the  immediate  area  change  at  the  throat  exit  produces  in  an  infinite  rate  of 
area  change  per  unit  length.  As  a  result,  this  abrupt  area  change  cannot  be  incorporated 
into  the  computational  domain.  To  simulate  the  flow  in  this  facility,  an  artificial  diverging 
nozzle  with  a  30°  half  angle  was  used  to  transition  from  the  outer  radius  of  the  throat  section 
to  the  outer  radius  of  the  expansion  section,  thereby  providing  a  finite  rate  of  area  change 
per  unit  length  throughout  the  computational  domain. 

The  second  significant  difference  between  the  2.44-m  blast  simulator  and  its  correspond¬ 
ing  computational  model  is  the  method  used  to  specify  the  initial  conditions.  As  stated 
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earlier,  the  high  pressure  driver  gas  and  the  ambient  air  are  separated  by  two  diaphragms 
that  have  an  intermediate  pressure  level  in  between  them.  In  the  computational  model,  the 
downstream  diaphragm  and  the  fluid  at  the  intermediate  pressure  state  are  completely  ig¬ 
nored.  Instead,  the  computational  domain  is  initialized  with  the  driver  gas  and  the  ambient 
air  immediately  adjacent  to  each  other  at  a  grid  point  corresponding  to  the  position  of  the 
upstream  diaphragm. 

The  model  of  the  2.44-m  blast  simulator  was  discretized  with  a  uniform  grid  spacing 
of  10  cm.  A  calculation  was  performed  to  simulate  a  test  that  had  been  conducted  in 
this  facility.  The  test  was  performed  with  an  ambient  pressure  of  101.5  kPa,  an  ambient 
temperature  of  313.7  K,  a  driver  overpressure  of  5.21  MPa,  and  driver  temperature  of  354  K. 
The  flow  histories  derived  from  the  experimental  data  are  compared  to  the  computed  histories 
in  Figures  6  through  8  and  summarized  in  Table  2. 


Table  2.  Summary  of  Results  for  2.44-m  Blast  Simulator  Test. 


Experiment 

Simulation 

Error 

Shock  Wave  Speed 

Incident  Static  Overpressure 

Dynamic  Pressure  Behind  Shock 

Mach  Number  Behind  Shock 

100.8  kPa 
31.3  kPa 
0.47 

480  m/s 
96.6  kPa 
28.9  kPa 
0.46 

-4.2% 

-7.7% 

-2.1% 

Static  Overpressure  Impulse  at  300  ms 
Dynamic  Pressure  Impulse  at  300  ms 

12.9  kPa-s 
2.08  kPa-s 

12.1  kPa-s 
2.39  kPa-s 

-6.2% 

-hl4.9% 

Figure  6.  2.44-m  Blast  Simulator  Test:  Static  Overpressure. 
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Mach  Number 


Figure  7.  2.44-m  Blast  Simulator  Test:  Dynamic  Pressure. 


Figure  8.  2.44-m  Blast  Simulator  Test:  Mach  Number. 


11 


The  static  overpressure  histories  are  provided  in  Figure  6  and  illustrate  the  experimental 
peak  static  overpressure  of  100.8  kPa.  The  shock  wave  speed  in  the  simulation  of  this 
experiment  was  480  m/s.  Using  the  Rankine-Hugoniot  relation  from  Equation  11,  the  peak 
static  overpressure  from  the  calculation  was  determined  to  be  96.6  kPa,  under-predicting 
the  experimental  peak  by  4.2%.  From  the  time  of  shock  arrival  to  about  350  ms,  the 
computational  static  overpressure  history  and  dynamic  pressure  history  of  Figure  7  exhibit  a 
“stair  step”  decay  pattern.  These  distinct  decreases  in  the  flow  histories  signify  the  passage 
of  expansion  waves  originating  from  the  driver  tube.  This  same  decay  pattern  can  be  seen 
in  the  experimental  data  in  the  two  figures  and  compare  well  with  the  computed  histories. 

In  shock  tube  tests  such  as  the  one  described  here,  when  the  primary  shock  reaches 
the  downstream  end  of  the  expansion  section,  it  encounters  the  sudden  area  change  at  the 
tunnel  exit.  No  longer  contained  by  the  walls  of  the  expansion  section,  the  shock  expands 
rapidly  into  the  surrounding  atmosphere.  The  sudden  expansion  of  the  shocked  gas  causes 
a  rarefaction  wave  to  form  at  the  tunnel  exit.  This  rarefaction  wave  travels  up  stream 
against  the  direction  of  fluid  flow.  The  pressure  on  the  downstream  side  of  the  rarefaction 
wave  is  lower  than  the  pressure  on  the  upstream  side.  As  the  rarefaction  wave  propagates  up 
stream,  the  fluid  particles  that  were  set  in  motion  by  the  primary  shock  encounter  the  reduced 
pressure  on  the  downstream  side  of  the  rarefaction  wave  and  are  accelerated.  The  code  is 
able  to  simulate  the  generation  of  the  rarefaction  wave  through  the  use  of  an  inflow-outflow 
boundary  condition  at  the  downstream  end  of  the  computational  domain.  This  boundary 
condition  forces  the  thermodynamic  pressure  at  the  outlet  to  be  always  equal  to  the  reference 
atmospheric  pressure.  In  Figure  6,  the  static  overpressure  decreases  suddenly  at  a  time 
of  about  350  ms  corresponding  to  the  arrival  of  the  rarefaction  wave  at  the  measurement 
position.  At  this  same  time,  the  dynamic  pressure  (see  Figure  7)  and  Mach  number  (see 
Figure  8)  are  observed  to  increase  as  a  result  of  the  acceleration  of  the  shocked  gas  by  the 
passage  of  the  rarefaction  wave.  The  use  of  the  inflow-outflow  boundary  condition  at  the 
downstream  end  of  the  computational  domain  allows  the  code  to  simulate  the  formation  and 
propagation  of  the  rarefaction  wave  as  evidenced  by  the  agreement  between  the  experimental 
and  computed  histories  after  the  arrival  of  the  rarefaction  wave  at  the  test  station. 

The  experimental  and  computational  static  overpressures  were  integrated  to  determine 
the  static  overpressure  impulse  and  the  diflferential  pressure  impulse.  At  300  ms,  the  static 
overpressure  impulse  of  the  experimental  data  has  reached  a  value  of  12.9  kPa-s,  while  the 
imp^ls0  from  the  computational  history  was  12.1  kPa-s,  a  diflPerence  of  6.2%.  In  the  case 
of  the  differential  pressure,  the  impulse  of  the  experimental  history  at  300  ms  is  2.13  kPa-s, 
while  that  of  the  simulation  is  2.48  kPa-s,  an  over-prediction  of  16.4%. 


6.3.  Variable  Area  Multiple  Driver  Blast  Simulator 

The  shock  tube  geometries  presented  until  now  have  been  relatively  simple  in  their 
configurations,  with  a  central  flow  regime  of  either  constant  or  variable  cross-sectional  area. 
The  last  level  of  geometric  complexity  to  be  studied  in  the  validation  of  the  code  is  that  in 
which  multiple  driver  tubes  feed  into  a  single  expansion  section.  Such  facilities  have  been 
used  in  Prance,  Germany,  and  the  United  States  to  simulate  exponentially  decaying  blast 
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flows  associated  with  the  detonation  of  tactical  nuclear  weapons.^  These  multi-driver  blast 
simulators  t3^ically  use  driver  tubes  of  assorted  lengths.  When  the  driver  tubes  are  charged 
to  a  uniform  initial  state  and  fired  simultaneously,  the  driver  tubes  produce  expansion  waves 
at  different  frequencies.  This  approach  helps  to  eliminate  the  “stair  step”  pattern  of  the 
pressure  histories,  allowing  them  to  be  more  like  an  exponential  decay. 

To  model  the  flow  in  multi-driver  shock  tubes  with  the  code,  the  actual  3-D  geometry  of 
such  a  facility  must  be  reduced  to  a  domain  with  one  computational  dimension  that  has  area 
changes  along  that  dimension.  This  is  accomplished  by  combining  the  total  available  flow 
area  into  an  equivalent  cross-sectional  area  at  every  point  along  the  length  of  the  facility. 
Using  this  technique,  the  1-D  computational  model  contains  the  same  total  volume,  mass, 
energy,  and  relative  area  changes  of  the  actual  facility  and  would  be  expected  to  produce 
adequate  results  for  longitudinal  positions  in  the  facility  where  the  flow  is  expected  to  be 
nearly  uniform  across  the  entire  cross  section.  However,  the  simplification  of  the  geometry 
into  this  1-D  form  makes  it  impossible  to  properly  study  the  flow  where  there  is  significant 
mixing  of  the  gas  exiting  the  different  driver  tubes. 

One  example  of  a  multi-driver  shock  tube  is  the  Large  Blast/Thermal  Simulator  (LB/TS), 
at  the  White  Sands  Missile  Range  in  New  Mexico,  and  operated  under  the  direction  of  the 
Defense  Special  Weapons  Agency  (DSWA).  A  photograph  of  this  facility  is  provided  in  Fig¬ 
ure  9.  The  LB/TS  has  nine  cylindrical  driver  tubes,  each  with  a  diameter  of  1.83  m.  These 
tubes  can  be  seen  near  the  lower  right  corner  of  the  photograph.  At  the  downstream  end  of 
each  driver  tube  is  a  converging  nozzle  1.45  m  in  length.  The  nozzle  transitions  the  avail¬ 
able  cross  section  from  the  driver  diameter  of  1.83  m  to  the  throat  diameter  to  0.91  m.  In 
the  throat  section,  the  distance  from  the  downstream  end  of  the  converging  nozzle  to  the 
diaphragm  is  1.63  m,  and  the  distance  from  the  diaphragm  to  the  downstream  end  of  the 
throat  section  is  2.29  m. 


The  lengths  of  each  of  the  cylindrical  driver  tubes  can  be  adjusted  by  securing  a  special 
plug  inside  each  tube  before  a  test.  With  these  plugs  placed  so  that  the  maximum  driver 
lengths  are  obtained,  the  lengths  of  the  cylindrical  driver  tubes  are  as  follow: 

•  two  tubes,  each  11.0  m  long, 

•  three  tubes,  each  19.9  m  long, 

•  two  tubes,  each  28.4  m  long,  and 

•  two  tubes,  each  36.2  m  long. 

The  expansion  section  is  170  m  long  and  is  in  the  shape  of  a  horizontal  semi-cylinder 
with  a  cross-sectional  area  of  163  m^.  The  upstream  end  of  the  expansmn  section  is  a 
vertical  wall,  with  the  downstream  ends  of  the  nine  throat  sections  flush  with  the  wall  At 
the  downstream  end  of  the  expansion  section  is  a  rarefaction  wave  eliminator  (RWE). 
The  RWE  can  make  the  expansion  tunnel  appear  to  be  infinitely  long  as  it  suppresses  the 
formation  of  an  expansion  wave  when  the  shocked  air  exits  the  expansion  section.  The  RWE 
is  mounted  on  railroad  tracks  so  that  it  can  be  moved  to  the  side  to  provide  equipment  access 
to  the  interior  of  the  expansion  section. 

A  calculation  was  performed  to  model  the  LB/TS  for  a  test  in  which  the  nine  driver 
tubes  had  been  configured  at  their  maximum  volume.  The  initial  overpressure  of  the  gas 
in  the  driver  tubes  was  3.44  MPa  and  the  initial  temperature  of  the  driver  gas  was  414  K. 
The  ambient  pressure  and  temperature  used  for  the  calculation  were  84.9  kPa  and  288.7  K, 
respectively.  The  geometry  of  the  upstream  end  of  the  computational  model  is  illustrated  in 
Figure  10.  The  model  was  developed  using  the  area  combination  method  described  earlier 
and  the  domain  was  discretized  with  a  uniform  grid  spacing  of  40  cm.  As  in  the  case  with  the 
2.44-m  blast  simulator,  an  artificial  diverging  nozzle  was  used  in  the  computational  model  to 
transition  the  area  from  the  exit  of  the  throat  sections  to  the  area  of  the  expansion  section. 
In  this  case,  a  diverging  nozzle  half-angle  of  45"  was  used. 


Figure  10.  LB/TS  Model  for  Single  Stage  Test:  Upstream  End. 

In  this  single  stage  experiment,  the  RWE  was  used  to  minimize  flow  disturbances  that 
would  be  created  by  the  sudden  expansion  of  the  shocked  gas  at  the  exit  plane  of  the 
expansion  section.  The  influence  of  the  RWE  on  the  exit  flow  is  replicated  in  the  calculation 
through  the  use  of  the  inflow/outflow  boundary  condition  described  earlier  and  changes  in 
the  grid  geometry  as  a  function  of  time.  A  converging  nozzle  is  placed  at  the  downstream 
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end  of  the  computational  domain  and  the  exit  area  is  changed  with  time  so  that  it  matches 
the  exit  area  history  of  the  RWE  in  the  experiment.  The  variation  in  grid  geometry  at  the 
downstream  end  is  illustrated  in  Figure  11.  This  figure  shows  the  different  exit  areas  at 
times  of  0.0,  0.6,  1.0,  and  1.4  s.  The  exit  area  is  held  constant  from  the  beginning  of  the 
calculation  until  the  arrival  of  the  shock  at  the  exit  plane  at  about  0.36  s.  After  then,  it  is 
decreased  during  the  event  to  accelerate  the  flow  and  reduce  the  effect  of  the  short  tunnel 
on  the  flow  at  the  test  section. 


Figure  11.  LB/TS  Model  for  Single  Stage  Test:  Downstream  End. 

The  LB/TS  is  configured  so  that  instrumentation  may  be  placed  at  a  number  of  positions 
throughout  the  expansion  section  for  measuring  blast  histories.  For  the  purpose  of  code 
validation,  comparisons  were  made  using  data  collected  at  a  position  which  is  105  m  from 
the  upstream  end  of  the  expansion  section.  This  corresponds  to  the  typical  location  of 
equipment  that  is  placed  in  the  expansion  tunnel  for  blast  testing.  The  computed  histories 
from  this  test  station  are  compared  to  the  experimental  data  in  Figures  12  through  14.  The 
computational  and  experimental  data  are  also  summarized  in  Table  3.  Figure  12  shows  a 
peak  static  overpressure  from  the  experiment  of  80.8  kPa  while  that  of  the  simulation  was 
75.9  kPa.  The  static  overpressure  positive  phase  ends  at  approximately  750  ms,  at  which 
time  the  impulse  was  20.4  kPa-s  for  the  calculation,  over-predicting  the  experimental  value 
of  18.8  kPa-s  by  8.5%. 

At  a  time  of  approximately  650  ms  in  Figure  12,  there  is  a  spike  in  the  computed  static 
overpressure  history,  followed  by  a  sudden  decrease  in  the  static  overpressure  of  both  the 
experimental  and  computed  histories.  At  this  same  time  in  the  dynamic  pressure  histo¬ 
ries  (see  Figure  13),  there  is  an  increase.  This  is  an  indication  that  the  open  area  at  the 
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static  Overpressure  (kPa) 


Figure  12.  LB/TS  Single  Stage  Test:  Static  Overpressure. 


Time  (s) 


Figure  13.  LB/TS  Single  Stage  Test:  Dynamic  Pressure. 


Figure  14.  LB/TS  Single  Stage  Test:  Mach  Number. 

downstream  end  of  the  expansion  tunnel  was  not  sufficiently  decreased  to  fully  eliminate  the 
formation  of  an  expansion  wave.  This  resulted  in  the  formation  of  a  weak  expansion  wave 
at  the  downstream  end  of  the  expansion  section  that  propagated  up  stream  and  caused  a 
slight  increase  in  the  particle  velocity,  resulting  in  the  pressure  drop  and  dynamic  pressure 
increase  evident  the  the  figures.  The  increase  in  particle  velocity  beginning  at  650  ms  is 
also  evident  in  the  Mach  number  history  of  Figure  14.  The  simulation  properly  predicts 
the  effect  of  the  expansion  wave  on  the  Mach  number  and  dynamic  pressure  histories  but 
over-predicts  its  magnitude.  After  the  arrival  of  the  expansion  wave,  the  simulated  Mach 
number  history  reaches  a  maximum  of  0.30  at  700  ms,  while  the  experimental  result  increases 
to  only  0.21.  At  the  same  times  the  simulated  dynamic  pressure  increases  to  6.35  kPa  while 
the  experiment  increases  to  3.38  kPa. 


Table  3.  Summary  of  Results  for  LB/TS  Single  Stage  Test. 


Experiment 

Simulation 

Error 

Shock  Wave  Speed 

Incident  Static  Overpressure 

Dynamic  Pressure  Behind  Shock 

Mach  Number  Behind  Shock 

80.3  kPa 
23.9  kPa 
0.45 

455  m/s 
77.3  kPa 
22.2  kPa 
0.44 

-3.7% 

-7.1% 

-2.2% 

Static  Overpressure  Impulse  at  750  ms 
Dynamic  Pressure  Impulse  at  1.0  s 

18.8  kPa-s 
4.16  kPa-s 

20.4  kPa-s 
4.81  kPa-s 

+8.5% 

+15.6% 
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Before  the  arrival  of  the  weak  expansion  wave  at  the  measurement  position,  the  computed 
dynamic  pressure  history  is  consistently  greater  than  the  measured  data  from  300  ms  to 
650  ms.  As  a  result  of  this  and  the  over-prediction  of  the  increased  particle  velocity  from  the 
weak  expansion  wave,  the  computed  dynamic  pressure  impulse  at  1.5  s  was  15.6%  greater 
than  that  obtained  from  the  experiment. 

7.  MULTIPLE  STAGE  SHOCK  TUBE  FLOW  SIMULATIONS 

Recently,  a  need  has  arisen  to  generate  complex  blast  flow  histories  in  shock  tubes. 
The  most  significant  characteristics  of  these  flow  histories  are  the  high  levels  of  dynamic 
pressure  and  the  presence  of  more  than  one  dynamic  pressure  peak  in  the  event.  For 
blast  simulators  with  multiple,  high  pressure  driver  tubes,  one  way  to  accomplish  this  is  to 
use  a  staged  opening  of  the  driver  tubes.  Using  this  technique,  a  shock  generated  by  the 
opening  of  a  set  of  drivers  can  be  used  to  reinforce  or  diminish  the  presence  of  a  particular 
flow  pattern  from  a  set  of  driver  tubes  opened  previously.  Flow  fields  of  this  type  can  be 
simulated  using  2-D  and  3-D  CFD  codes,  but  the  computational  analysis  that  is  required 
can  be  very  time  consuming  as  the  simulation  of  one  multi-stage  test  may  actually  consist 
of  several  calculations  that  provide  information  to  each  other.  To  decrease  the  time  spent  in 
designing  a  multi-stage  shock  tube  test,  the  capability  to  model  this  type  of  event  was  added 
to  the  code.  This  1-D  modeling  capability  allows  the  analyst  to  study  a  range  of  possible  test 
conditions  in  a  short  period  of  time.  This  analysis  would  eliminate  the  cases  that  have  little 
chance  of  success  and  allow  the  analyst  to  perform  more  detailed  2-D  or  3-D  computations 
for  those  cases  that  have  a  greater  chance  of  success. 

Because  the  code  solves  the  equations  of  motion  in  only  one  dimension,  it  is  not  possible 
for  a  single  computational  model  to  replicate  the  individual  driver  tubes  firing  at  diflferent 
times.  To  make  the  simulation  possible,  separate  calculations  are  performed  to  model  the 
exit  flow  from  the  banks  of  driver  tubes  for  each  individual  stage  of  a  multi-stage  event.  For 
each  of  these  independent  calculations,  the  primitive  flow  parameters  (pressure,  density,  and 
velocity)  are  collected  as  a  function  of  time  at  the  driver  exit.  The  driver  exit  history  data 
are  then  processed  and  used  as  an  inlet  boundary  condition  for  a  separate  simulation  of  the 
multi-stage  event. 

Since  only  one  computational  dimension  exists,  the  multiple  driver  exit  histories  may 
not  be  used  independently  to  drive  the  simulation  of  the  multi-stage  flow.  The  histories  must 
be  combined  and  provided  as  a  single  flow  history.  Conservation  laws  are  used  to  provide  a 
basis  for  combining  the  multiple  histories  into  a  single  history  of  primitive  flow  variables  and 
inlet  area.  Specifically,  the  total  mass,  momentum,  energy,  and  inlet  cross-sectional  area  are 
conserved  in  the  combined  history  as  defined  in  Equations  15  and  16. 

i^i)net  =  E  (9i)m 
m=l 

n 

^net  ~ 

m=l 
18 


(16) 


In  these  equations,  qi  represents  the  conservation  parameter  of  mass,  momentum,  or  energy, 
and  A  represents  the  driver  exit  area.  The  summation  of  m  from  1  to  n  is  used  to  denote 
the  number  of  active  inlet  histories  at  the  point  in  time  when  the  calculation  is  performed. 
The  summations  produce  the  combined  inlet  flow  history  denoted  with  the  net  subscript. 

After  the  total  mass,  momentum,  energy,  and  inlet  cross-sectional  area  are  computed,  the 
net  density  is  computed  by  dividing  the  total  area  into  the  net  mass  term  as  in  Equation  17: 
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The  net  velocity  is  computed  by  dividing  the  mass  term  into  the  momentum  term  as  illus¬ 


trated  in  Equation  18: 
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From  Equations  3,  5,  and  15,  the  total  energy  from  the  individual  histories  is  given  by 
Equation  19: 


Using  the  net  area,  density  and  velocity  from  the  previous  equations.  Equation  19  can  be 
solved  for  the  net  pressure  as  shown  in  Equation  20: 
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Performing  this  set  of  calculations  produces  the  history  of  net  pressure,  density,  and 
velocity  which  is  used  as  a  time-dependent,  inlet  boundary  condition  for  the  multi-staged 
simulation. 


7.1.  Two-Stage  Tests 

Generation  of  high  dynamic  pressure  flows  with  a  two-stage  driver  firing  sequence  is 
accomplished  by  first  releasing  the  gas  from  a  specified  bank  of  drivers.  When  the  leading 
shock  from  this  first  bank  reaches  the  downstream  end  of  the  expansion  section,  a  rarefaction 
wave  is  formed,  which  travels  up  stream  in  the  expansion  section,  accelerating  the  shocked 
gas  as  it  propagates  up  stream.  The  timing  of  flow  initiation  from  the  second  driver  bank 
is  such  that  the  shock  formed  by  the  second  driver  release  arrives  at  the  desired  location  at 
the  same  time  as  the  rarefaction  wave.  This  coincidental  arrival  of  the  second  shock  and  the 
rarefaction  wave  produces  a  sudden  increase  in  the  dynamic  pressure. 
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Two  tests  of  this  type  have  been  performed  in  the  LB/TS.  The  RWE  was  not  used  in 
either  test.  In  both  cases,  the  first  bank  of  drivers  was  comprised  of  the  five  longest  individual 
driver  tubes,  having  a  combined  total  volume  385  m^.  The  second  bank  was  comprised  of  the 
four  shortest  driver  tubes  with  a  total  volume  of  168  m®.  Releasing  the  large  bank  first  forces 
the  majority  of  the  driver  gas  into  the  expansion  tunnel  early  in  the  event.  The  purpose  of 
the  second  bank  is  to  propagate  a  shock  behind  the  driver  gas  released  by  the  first  bank  and 
drive  it  down  the  tunnel. 

Simulation  of  a  two-stage  test  is  accomplished  through  the  setup  and  execution  of  three 
separate  calculations.  The  first  calculation  simulates  the  independent  firing  of  the  first  bank 
of  driver  tubes  into  the  expansion  section.  The  same  is  done  for  the  second  bank  in  a 
separate  calculation.  In  each  of  these  two  calculations,  the  fiow  histories  are  recorded  at  the 
downstream  end  of  the  throat  section,  where  the  driver  gas  flows  into  the  expansion  section. 
These  two,  independent  flow  histories  are  then  combined  using  the  methodology  described 
earlier  in  this  section  to  produce  a  single  representative  history  of  pressure,  density,  velocity, 
and  inlet  area. 

The  final  calculation  is  designed  to  simulate  only  that  portion  of  the  shock  tube  down 
stream  from  the  throat  section.  At  the  upstream  end  of  this  computational  model,  the  area 
is  varied  with  time  to  simulate  the  increase  in  available  flow  area  when  the  second  driver 
bank  is  released.  The  combined  history  from  the  individual  driver  calculations  is  supplied 
as  a  forcing  function  at  the  upstream  boundary. 

In  the  first  of  these  two-stage  tests,  the  bank  of  five  drivers  was  charged  to  an  overpressure 
of  2.96  MPa  and  a  temperature  of  310.9  K,  while  the  bank  of  four  drivers  was  initialized 
with  an  overpressure  and  temperature  of  6.41  MPa  and  449.8  K,  respectively.  When  the  test 
was  performed,  one  of  the  diaphragms  in  the  five-driver  bank  failed  to  open,  leaving  only 
four  tubes  to  drive  the  first  pulse.  The  second  driver  bank  was  fired  300  ms  after  the  first. 

To  simulate  this  test,  the  computational  model  of  the  first  driver  bank  consisted  only  of 
the  four  driver  tubes  that  opened.  The  model  of  the  second  driver  bank  was  also  developed 
using  the  area  combination  method  described  earlier.  These  two  calculations  were  performed 
and  the  flow  histories  at  the  throat  exit  were  used  to  generate  the  inlet  history  for  the  staged 
flow  calculation.  In  the  staged  flow  calculation,  the  inlet  area  was  varied  with  time  to  model 
the  changing  available  flow  area  as  the  drivers  are  opened.  Prom  a  time  of  zero  to  300  ms,  the 
inlet  cross-sectional  area  was  2.63  m^,  corresponding  to  the  area  of  the  four  active  drivers 
in  the  first  bank.  At  300  ms,  the  inlet  area  was  immediately  doubled  to  account  for  the 
additional  four  drivers  in  the  second  bank,  then  held  constant  for  the  remainder  of  the 
simulation.  Figure  15  shows  the  pressure,  density,  and  velocity  histories  of  the  inlet  flow.  In 
these  figures,  a  time  of  zero  corresponds  to  the  time  of  diaphragm  rupture  of  the  first  driver 
bank.  Immediately  after  the  diaphragm  rupture  of  the  first  bank,  the  pressure,  density,  and 
velocity  increase  and  then  begin  to  decay  gradually.  At  300  ms,  the  inlet  area  is  immediately 
opened  to  the  equivalent  cross-sectional  area  of  all  eight  throat  sections,  causing  the  sharp 
drop  in  the  inlet  histories.  Immediately  after  this  increase  in  inlet  area,  the  flow  from  the 
second  driver  bank  is  initialized  and  the  fiow  parameters  again  increase  suddenly,  then  slowly 
decay  as  the  driver  tubes  empty. 
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Figure  15.  LB/TS  Two-Stage  Test  #1:  Inlet  Histories. 
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A  static  overpressure  history  was  recorded  at  a  position  103.8  m  from  the  upstream  end 
of  the  expansion  tunnel,  stagnation  pressure  measurements  were  recorded  at  105.8  m  and 
141.1  m,  and  differential  pressure  data  were  collected  at  141.2  m.  The  conditions  behind 
the  first  shock  and  the  impulse  values  for  the  experimental  and  computed  histories  are 
provided  in  Table  4.  The  static  overpressure  data  are  plotted  in  Figure  16  and  the  stagnation 
overpressure  history  data  from  the  gauge  at  105.8  m  are  shown  in  Figure  17.  These  figures 
show  the  arrival  of  the  first  shock  at  260  ms,  and  the  arrival  of  the  second  shock  at  520  ms. 
The  calculation  under- predicts  the  amplitude  of  both  the  incident  static  overpressure  and  the 
incident  stagnation  pressure.  The  arrival  of  the  rarefaction  wave  at  the  103.8-m  location  is 
visible  in  both  the  experimental  and  computed  static  overpressure  histories  at  approximately 
700  ms,  where  the  downward  slopes  of  the  curves  increase  and  lead  to  the  the  end  of  the 
positive  phase  of  the  event  at  approximately  800  ms. 


Table  4.  Summary  of  Results  for  LB/TS  Two-Stage  Test 


Experiment 

Simulation 

Error 

Shock  Wave  Speed 

Incident  Static  Overpressure 

Dynamic  Pressure  Behind  Shock 

Mach  Number  Behind  Shock 

40.9  kPa 
6.6  kPa 
0.27 

395  m/s 
34.2  kPa 
4.7  kPa 
0.24 

-16.4% 

-28.8% 

-11.1% 

Static  Overpressure  Impulse  at  800  ms  (103.8  m) 
Static  Overpressure  Impulse  at  900  ms  (141  m) 
Dynamic  Pressure  Impulse  at  1.3  s  (141  m) 

17.7  kPa-s 
11.0  kPa-s 
8.0  kPa-s 

17.0  kPa-s 
9.0  kPa-s 
8.7  kPa-s 

-4.0% 

-18.2% 

+8.8% 

station  Location  103.8  m 


Figure  16.  LB/TS  Two-Stage  Test  #1:  Static  Overpressure  at  103.8  m. 
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The  stagnation  overpressure  history  recorded  at  the  141. 1-m  location  and  the  differen¬ 
tial  pressure  history  recorded  at  141.2  m  were  combined  to  produce  the  static  overpressure, 
dynamic  pressure,  and  Mach  number  histories  in  Figures  18  through  20.  These  measure¬ 
ments  were  recorded  at  a  location  farther  down  stream  in  the  tunnel  than  those  provided  in 
Figures  16  and  17  and  show  that  the  rarefaction  wave  arrives  at  the  gauge  location  shortly 
before  the  second  shock.  The  arrival  of  the  rarefaction  wave  is  marked  by  the  sudden  de¬ 
crease  in  the  static  overpressure  history  at  530  ms  and  the  increase  in  the  dynamic  pressure 
and  Mach  number  histories  at  the  same  time.  The  second  shock  then  arrives  at  600  ms  and 
combines  with  the  effect  of  the  rarefaction  wave  to  create  a  second  peak  in  the  differential 
pressure  history  of  approximately  72  kPa.  The  positive  phase  of  the  static  overpressure  ends 
at  approximately  825  ms  for  both  the  experiment  and  the  calculation.  In  the  negative  phase 
of  the  static  overpressure  history,  a  minimum  pressure  is  reached  at  approximately  975  ms, 
after  which  time,  the  histories  asymptotically  approach  zero. 

In  the  dynamic  pressure  and  Mach  number  histories  in  Figures  19  and  20,  the  compu¬ 
tational  results  closely  follow  the  experimental  data  after  the  arrival  of  the  second  shock 
until  approximately  800  ms.  Between  800  ms  and  950  ms,  the  computed  dynamic  pressure 
reaches  a  value  of  15  kPa  while  the  experimental  data  during  this  same  period  is  approx¬ 
imately  11  kPa.  The  Mach  number  shows  a  similar  over-prediction  during  this  period  of 
time. 

The  second  two-stage  test  performed  in  the  LB/TS  was  performed  with  all  the  drivers 
charged  to  an  overpressure  of  6.2  MPa  and  a  temperature  of  310.9  K.  Like  the  previous 
test  discussed,  this  test  was  run  with  the  bank  of  the  five  largest  drivers  fired  first,  followed 
300  ms  later  by  the  bank  of  the  four  smallest  drivers.  In  this  case,  all  the  drivers  opened 
properly.  The  experiment  was  performed  with  a  static  overpressure  gauge  at  the  103. 8-m 
location,  stagnation  overpressure  gauges  at  105.8  m  and  155.6  m,  and  differential  pressure 
gauges  at  105.8  m,  141.2  m,  and  144.5  m. 

The  computed  and  simulated  static  overpressures  for  the  103.8-m  location  are  compared 
in  Figure  21.  The  incident  static  overpressure  produced  by  the  simulation  was  28.2%  less  than 
that  of  the  experiment.  This  lower  shock  amplitude  produced  by  the  simulation  resulted  in  a 
difference  in  shock  arrival  time  of  approximately  8  ms.  Despite  this  late  arrival  of  the  incident 
shock,  the  simulation  properly  predicts  the  arrival  time  and  amplitude  of  the  second  shock. 
From  the  arrival  of  the  second  shock  until  the  end  of  the  positive  phase  at  approximately 
825  ms,  the  computational  result  closely  follows  the  experimental  data. 

The  experimental  histories  from  the  two  gauges  at  the  105. 8-m  location  were  used  to 
compute  the  static  overpressure,  dynamic  pressure,  and  Mach  number  histories  that  are 
provided  in  Figures  22,  23,  and  24,  respectively.  Like  the  history  at  the  103.8-m  location, 
these  figures  show  the  under-prediction  of  the  amplitude  of  the  incident  shock,  the  late 
arrival  of  the  computed  incident  shock,  and  the  correct  prediction  of  the  arrival  time  of  the 
second  shock.  After  the  arrival  of  the  second  shock,  the  dynamic  pressure  and  Mach  number 
histories  pass  through  the  scatter  in  the  experimental  data.  There  is  a  sudden  decrease  in  the 
experimental  dynamic  pressure  and  Mach  number  at  approximately  1.03  s.  This  decrease  is 
reproduced  in  the  simulation  but  occurs  approximately  10  ms  earlier. 
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re  23.  LB/TS  Two-Stage  Test  #2:  Dynamic  Pressure  at  105.8  m. 


Station  Location  105.8  m 


Time  (s) 


Figure  24.  LB/TS  Two-Stage  Test  #2:  Mach  Number  at  105.8  m. 
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25.  LB/TS  Two-Stage  Test  #2:  Differential  Pressure  at  141.2  m 


Station  Location  1 44.5  m 


Time  (s) 


Figure  26.  LB/TS  Two-Stage  Test  #2:  Differential  Pressure  at  144.5  m 


station  Location  155.6  m 


Figure  27.  LB/TS  Two-Stage  Test  #2;  Stagnation  Overpressure  at  155.6  m. 


Table  5.  Summary  of  Results  for  LB/TS  Two-Stage  Test  #2. 


Experiment 

Simulation 

Error 

Shock  Wave  Speed 

- 

424  m/s 

- 

Incident  Static  Overpressure 

75.3  kPa 

54.1  kPa 

-28.2% 

Dynamic  Pressure  Behind  Shock 

21.2  kPa 

11.3  kPa 

-46.7% 

Mach  Number  Behind  Shock 

0.43 

0.34 

-20.9% 

Static  Overpressure  Impulse  at  800  ms  (103.8  m) 

24.9  kPa-s 

23.6  kPa-s 

-5.2% 

Static  Overpressure  Impulse  at  800  ms  (105.8  m) 

27.6  kPa-s 

23.0  kPa-s 

-16.7% 

Dynamic  Pressure  Impulse  at  1.4  s  (105.8  m) 

12.7  kPa-s 

12.0  kPa-s 

-5.5% 

Differential  Pressure  Impulse  at  1.4  s  (141.2  m) 

14.5  kPa-s 

15.2  kPa-s 

-f-4.8% 

Differential  Pressure  Impulse  at  1.4  s  (144.5  m) 

15.3  kPa-s 

15.4  kPa-s 

-hO.7% 

Stagnation  Overpressure  Impulse  at  1.4  s  (155.6  m) 

25.7  kPa-s 

19.3  kPa-s 

-24.9% 
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The  remaining  differential  pressure  and  stagnation  overpressure  histories  are  shown  in 
Figures  25  through  27.  Because  these  histories  were  collected  at  a  location  farther  down 
stream,  the  difference  in  shock  arrival  time  between  the  experiment  and  the  calculation  is 
greater  than  it  was  at  the  105. 8-m  location,  22  ms  to  8  ms.  Similar  to  the  first  two-stage  test, 
the  histories  collected  close  the  exit  plane  of  the  expansion  section  show  the  rarefaction  wave 
arriving  at  the  gauge  locations  shortly  before  the  second  shock.  Like  the  results  from  the 
upstream  locations,  the  computation  under-predicts  the  flow  conditions  behind  the  incident 
shock  but  properly  predicts  the  arrival  time  of  the  second  shock  and  the  flow  histories  that 
follow  the  second  shock. 


7.2.  Four-Stage  Test 

The  primary  purpose  of  the  two-stage  staggered  driver  testing  that  has  been  presented 
here  is  to  produce  blast  histories  with  high  levels  of  dynamic  pressure.  Increasing  the  initial 
driver  pressure  results  in  increased  peak  dynamic  pressure  and  dynamic  pressure  impulse. 
However,  this  increase  in  driver  pressure  also  results  in  a  leading  shock  of  high  amplitude. 
When  testing  military  equipment  in  such  a  high  dynamic  pressure  environment,  a  strong 
leading  shock  may  cause  unwanted  structural  damage  to  the  target.  To  decrease  the  likeli¬ 
hood  of  causing  structural  damage  with  a  strong  leading  shock,  it  was  postulated  that  the 
driver  gas  could  be  released  with  a  larger  number  of  driver  flow  stages  in  the  event.  This 
would  have  the  effect  of  reducing  the  strength  of  each  individual  shock  and  diffusing  the 
initial  shock  energy  delivery  over  a  set  of  incremental  pressure  increases  rather  than  a  single, 
strong  shock. 

Such  a  test  was  performed  in  the  LB/TS.  In  this  test,  the  primary  measurements  were 
made  at  positions  of  103.8  m  and  105.8  m  from  the  upstream  end  of  the  expansion  tunnel. 
The  test  consisted  of  four  banks  of  drivers  in  which  the  initial  overpressure  and  temperature 
were  5.2  MPa  and  288.7  K,  respectively.  Like  the  two-stage  tests,  the  RWE  was  not  used 
in  order  to  take  advantage  of  the  rarefaction  wave  that  would  accelerate  the  flow  in  the 
expansion  tunnel  and  increase  the  dynamic  pressure.  The  flow  field  was  created  by  first 
firing  the  two  36.2-m  driver  tubes  (t=0  s),  then  the  two  11.0-m  tubes  at  t=100  ms,  followed 
by  the  three  19.9-m  tubes  at  t=225  ms,  and  finally  the  two  28.4-m  tubes  at  t=425  ms.  With 
this  configuration,  the  first  two  stages  of  flow  were  intended  to  provide  shaping  of  the  early 
time  history  and  introduce  dense  driver  gas  into  the  expansion  tunnel.  Then,  the  last  two 
stages  were  used  to  further  accelerate  the  gas  released  by  the  first  two  stages. 

Just  as  the  two-stage  experiment  required  three  calculations  to  simulate  the  event,  this 
four-stage  test  required  a  total  of  five  calculations  to  produce  the  overall  simulation.  Four 
independent  calculations  were  performed  to  simulate  the  flow  from  each  of  the  driver  banks 
emptying  into  the  expansion  section.  The  flow  histories  from  these  four  throat  sections  were 
collected  and  combined  using  the  method  described  earlier  to  generate  the  inlet  history  in 
Figure  28. 

These  histories  were  then  used  as  a  forcing  function  at  the  inlet  boundary.  The  cross- 
sectional  area  of  the  inlet  boundary  of  the  computational  domain  was  varied  with  time  to 
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Figure  28.  LB/TS  Four-Stage  Test:  Inlet  Fluid  Histories. 
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replicate  the  changes  in  area  associated  with  the  increasing  number  of  drivers  during  the 
event.  In  Figure  29,  the  inlet  boundary  area  is  plotted  as  a  function  of  time  to  show  the 
variation  in  the  geometry  during  the  calculation.  During  the  first  100  ms,  the  inlet  cross- 
sectional  area  was  1.31  m^,  the  sum  of  the  cross-sectional  areas  of  the  two  throat  sections. 
The  area  is  then  increased  at  the  appropriate  times  to  provide  the  equivalent  cross-sectional 
area  of  four,  seven,  and  finally  nine  throat  sections.  After  the  last  bank  of  drivers  begins  to 
fire,  the  inlet  cross-sectional  area  is  held  constant  for  the  remainder  of  the  calculation. 
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Figure  29.  LB/TS  Four-Stage  Test:  Inlet  Area  History. 

A  static  overpressure  gauge  was  used  at  the  103. 8-m  location  and  its  recorded  history 
is  compared  to  the  simulation  in  Figure  30.  In  this  figure,  four  distinct  increases  in  the 
experimental  and  computed  histories  signify  the  arrival  of  the  shocks  from  the  individual 
driver  pulses.  The  computation  under-predicts  the  incident  static  overpressure  by  22%  and 
is  the  cause  of  the  difference  in  arrival  times  of  the  first  shock.  The  arrival  times  of  the 
second,  third,  and  fourth  shocks  in  the  simulation  match  those  of  the  experiment.  The  static 
overpressure  positive  phase  for  the  experiment  ends  at  approximately  840  ms,  while  that  of 
the  calculation  ends  at  about  870  ms. 

A  stagnation  overpressure  gauge  and  a  differential  pressure  gauge  were  employed  at  the 
105.8-m  location  and  were  used  to  generate  the  static  overpressure,  dynamic  pressure,  and 
Mach  number  histories  of  Figures  31,  32,  and  33.  As  in  the  case  of  the  static  overpressure, 
the  computed  histories  at  the  105.8-m  location  exhibit  the  dominant  flow  characteristics 
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Figure  30.  LB/TS  Four-Stage  Test:  Static  Overpressure  at  103.8  m 


Figure  31.  LB/TS  Four-Stage  Test:  Static  Overpressure  at  105.8  m. 


observed  in  the  experimental  data.  Table  6  provides  additional  summary  information  from 
the  experimental  and  computational  histories. 


Table  6.  Summary  of  Results  for  LB/TS  Four-Stage  Test. 


Experiment 

Simulation 

Error 

Shock  Wave  Speed 

Incident  Static  Overpressure 

Dynamic  Pressure  Behind  Shock 

Mach  Number  Behind  Shock 

40.0  kPa 
6.3  kPa 
0.27 

391  m/s 
31.2  kPa 
3.9  kPa 
0.22 

-22.0% 

-38.1% 

-18.5% 

Static  Overpressure  Impulse  at  850  ms  (103.8  m) 
Static  Overpressure  Impulse  at  850  ms  (105.8  m) 
Dynamic  Pressure  Impulse  at  1.4  s  (105.8  m) 

20.2  kPa-s 
22.0  kPa-s 
8.9  kPa-s 

20.5  kPa-s 
20.0  kPa-s 
10.4  kPa-s 

+1.5% 

-9.1% 

-bl6.9% 

8.  SUMMARY 

This  report  has  documented  the  development  and  validation  of  a  computational  tool  for 
simulating  complex  flow  flelds  in  shock  tubes  and  blast  simulators.  The  tool  is  easy  to  use 
and  requires  only  20  to  30  lines  of  input  data  to  completely  define  a  geometry,  boundary 
conditions,  initial  fluid  conditions,  and  run  control  information.  Most  of  the  calculations 
illustrated  in  this  report  required  only  a  few  minutes  of  run  time  to  complete  on  an  average 
workstation.  The  most  complex  staggered  flow  calculation  ran  in  about  1  hour.  With  its 
ease  of  use  and  quick  turn-around  time,  this  code  makes  it  possible  to  perform  complex  flow 
studies  consisting  of  numerous  calculations  in  a  short  period  of  time.  These  traits  make 
the  code  well  suited  as  a  design  tool  for  shock  tube  experiments  or  for  deriving  first  order 
estimates  of  flow  fields  for  large  3-D  computational  studies. 

The  variety  of  shock  tube  geometries  presented  here  provides  an  indication  of  the  capa¬ 
bilities  and  limitations  of  the  code.  Overall,  the  code  performed  best  in  simulating  shock 
tubes  with  simple  geometries.  Differences  between  computed  and  experimental  flow  histories 
were  greatest  when  simulating  the  most  complex  geometries.  For  the  single  stage  flow  cases, 
the  predictions  of  the  incident  static  overpressure  were  typically  within  10%  of  the  experi¬ 
mental  results.  For  these  same  cases,  the  difference  between  the  simulated  and  experimental 
static  and  dynamic  impulses  was  less  than  20%. 

For  multi-staged  shock  tube  flows,  the  code  inherently  under-prediced  the  amplitude 
of  the  leading  shock  but  generally  followed  the  dominant  trends  in  the  experimental  flow 
histories.  The  error  in  the  prediction  of  the  incident  static  overpressure  was  greatest  for 
the  second  two-stage  test  and  was  28.2%.  The  difference  in  impulse  values  between  the 
experiment  and  computation  was  typically  less  than  15%  for  the  multi-staged  flows.  The 
code  demonstrated  that  it  can  accurately  reproduce  the  arrival  times  of  multiple  shocks  in 
a  multi-stage  event,  making  it  a  useful  tool  for  designing  experiments  that  require  such  flow 
characteristics. 
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